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Summary. Propagation of monochromatic elastic waves across cracks is investigated in ID, 
both theoretically and numerically. Cracks are modeled by nonlinear jump conditions. The 
mean dilatation of a single crack and the generation of harmonics are estimated by a pertur- 
bation analysis, and computed by the harmonic balance method. With a periodic and finite 
network of cracks, direct numerical simulations are performed and compared with Bloch- 
Floquet's analysis. 



1 Introduction 

Failure processes resulting in a crack generally produce rough crack faces. Once crack open- 
ing has taken place, and the crack faces have undergone slight relative sliding displacement, 
the crack will never completely close again due to the nonconforming surfaces in partial con- 
tact. A complicated interaction between crack faces is expected, depending strongly on the 
magnitudes of the tractions transmitted across the rough surfaces in contact [9]. 

The interaction of ultrasonic waves with cracks has been investigated by many authors, 
assuming that the wavelength is much larger than a characteristic length of the roughness of 
the contacting surfaces. Linear slip-displacement models of crack-face interaction have been 
widely used [12, 10]. However, a non-physical penetration of contacting surfaces may occur in 
linear models. Moreover, laboratory experiments have shown that methods of non-destructive 
evaluation based on linear models may fail to detect partially closed cracks [13]. 

Here, we study wave propagation with a nonlinear model of contact proposed in [1, 2]. A 
monochromatic compressional wave propagates normally to a plane flaw surface, leading to a 
ID problem detailed in section 2. Analysis of scattered fields is performed in section 3. With 
a single crack, the generation of harmonics and the mean dilatation of the crack are addressed 
analytically and numerically. Propagation through periodic networks of contact nonlinearities 
is studied by Bloch-Floquet's analysis [8] and simulations. Numerical experiments are pro- 
posed in section 4. Conclusions are drawn and future perspectives are suggested in section 
5. 
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2 Problem statement 
2.1 Configuration 




Fig. 1. Elastic media Qq and Qi separated by vacuum, with rough contact surfaces. Static 
(left) and dynamic (right) case, with incident (I), reflected (R) and transmitted (T) waves. 



We consider a single crack with rough faces separating two media £2q and fl± linearly 
elastic and isotropic, with density p and elastic speed of the compressional waves c. These 
parameters are piecewise constant and may be discontinuous around the crack: (po, Co) if x € 
Qq, (pi , c\ ) if x € il\ . The media are subject to a constant static stress p. At rest, the distance 
between planes of average height is £,o(p) > (figure 1, left). An incident monochromatic 
wave, emitted by a ponctual stress source at x = x s in Qq, gives rise to reflected (in Qq) and 
transmitted (in Q t ) compressional waves. These perturbations in Qq and il\ are described by 
the ID elastodynamic equations 

dv da da 2 dv „ 2 c-/ % • 

Pjj = ~fa' -^j = pc 2 po c v 5 (x-x s ) sin cot, (1) 

where vo is the amplitude of the incident elastic velocity, and ft) = 2nf is the angular fre- 
quency of the source. The elastic velocity v = ^7, the elastic displacement u, and the elastic 
stress perturbation o" around p, are averaged fields per unit area in crack's plane. The dynamic 
stresses induced by the elastic waves affect the thickness ^(t) of the crack (figure 1, right). 
The constraint 

£=£o + M>£o-S>0 (2) 

must be satisfied, where [u] = u + — iT is the difference between the elastic displacements 
on the two sides of the crack, and 8(p) > is the maximum allowable closure [2], We also 
assume that the wavelengths are much larger than t, , neglecting the propagation time through 
the crack, and replacing it by a zero-thickness interface at x = a: [u] = [u(a, t)} = u{a + , t) — 
u(a~,t). 
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Fig. 2. Nonlinear relation between stress and slip displacement. 



2.2 Contact law 

Cracks are classically modeled by linear jump conditions [12] with stiffness K: 



[o(a,t)] = 0, 



[u(a, t)] = — a(a ± , t) 
K 



(3) 



If K — > +°°, welded conditions are recovered. Conditions (3) violate (2) under large compres- 
sion loadings: o(a , t) < —K8 => | < |q — 5. Hence, the linear regime induced by (3) is 
realistic only with very small perturbations. With larger ones, nonlinear jump conditions are 
required. 

Compression loading increases the number and the surface of contacting faces. Conse- 
quently, a smaller stress is needed to open than to close a crack of a given displacement; 
an infinite stress is even required to close the crack faces completely. This behavior can be 
modeled by the global jump conditions proposed in [1, 2] 



[o(a,t) 



= 0, 



[u(a,t) 



1 



a(a ± ,f) 



K l- ( j(a ± ,t)/(K8)' 



(4) 



satisfying (2) and implying a(a , t) <K8. The second relation in (4) is sketched in figure 2. 
The straight line with a slope K tangential to the hyperbola at the origin amounts to the linear 
jump conditions (3); as deduced from (4), the linear regime is valid only if ^(a*, t) \ < KS. 
A second limit-case not investigated here is obtained if 8 — > and K bounded: the hyperbola 
tends towards the nondifferentiable graph of the unilateral contact, denoted by bold straight 
segments in figure 2. 



3 Analysis of scattered fields 
3.1 Single crack 

Analytical approach. The scattered fields can be expressed in terms of Y(t) = [u(a, t)]. Fol- 
lowing [11,4] gives the nonlinear ordinary differential equation (ODE) 
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dY Y 

■d7 +p TTYj8= 2vQSmcot > (5) 

satisfied by Y, with j5 = K((pqcq)~ 1 + (pi cj) -1 ). Inspection of (5) and dimensional analysis 
show that Y is 5 times a function of |r and cot. To solve (5), we assume \Y\/8 < 1, which 
leads to the series 

~T7 + P ~&f~ ¥ " +1 = 2v <> sinfflr - (6) 
An approximate solution of (6) is sought by a perturbation approximation (PA) [3] 

Y(t) = £ Y„(t) = £ JLy B (0. (7) 
n=0 n=0 ° 

Nothing ensures that (7) converges: |y|/5 < 1 is always satisfied when Y < 0, but it is true 
when y > only if ^ is sufficiently small. Plugging (7) into (6) and identifying the terms 
with identical power of 8 leads to an infinite series 

-jj- +j8>'0 = 2v sin cm, 

^+/3>'i=/3>i (8) 

dyi 
dt 



^+j8>'2 = /3>-o(2yi->-g), 



and so on. There is no influence of the truncation order N on the accuracy of the solution Y n 
with n < N. The recursive and linear ODE are much simpler to solve than the nonlinear ODE 
(5), even if computing Y n with n > 2 is cumbersome. This computation has been automatized 
with computer algebra tools. Since the number of terms in Y„ is roughly 2 2n+1 , very high 
orders are currently out of reach. Computing Y n up to N = 8 takes 20 mn on a Pentium 3 GHz. 
We detail the case N = 1. Setting <f>o = arctan j, the zero-th order periodic solution of (7)-(8) 
is 

y (f)=2^— 1 sin(fl»-m). (9) 

P v/l + (ffl/j8) 2 

Setting <f>] = arctan ^j, the first-order periodic solution of (7)-(8) is 

y lW=2^ / 2 fl- . 1 =sin(2fl»-2m + yi)). (10) 
/3 2 5 l + (co//3) 2 ^ ^ 1 + (2ffl //3) 2 / 

Two properties are deduced from (10). First, the mean value of Y\ is non-null 

Fi=2-^- i y>0. (11) 

(8 2 5 l + (ffl//3) 2 

_ / , \2n+2 _ 

More generally, the mean value y2«+i of y2„+i is proportional to 8 ( J , and Yin = 0. 

The mean thickness of the crack £ deduced from (2) satisfies t, = t,Q + Y > £0, with 
y = £y„. The dilatation predicted here is similar to the DC signal measured by [6] with 
glass-piezoceramic interface. The second property deduced from (10) concerns the term with 
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angular frequency 2 a. Its importance is quantified by the ratio of amplitudes between sinu- 
soidal terms in (10) and (9) 

n = ^ , \ (12) 

p d ^i+( £0 //3)yi+( 2ffl / J 8) 2 

Consequently, nonlinear effects increase with and decrease with j . 

Numerical approach. To compute the scattered fields with high accuracy and no limi- 
tation about the range of validity, one implements the numerical harmonic balance method 
(HBM). The periodic elastic displacements are written {ko = G)/cq, k\ = (0/c\) 

u;(x, t) = — {cos(ft>r — kox) — 1} , 

u R (x, t) =Rq+ ^ ^R"smn(at + k x)+R b 1 cosn(at + k x)^, ^) 
u T (x,t) = Tq + £ sinn(at-k\x) + T b cosn (at - k\ x)| . 

n— 1 

The elastic stresses are deduced from (13). Fields are truncated at N and injected in (4). The 
truncation implies that Rn' b and T"' b depend on N. The terms with identical trigonometric 
arguments are put together, and the terms with a trigonometric argument greater than N at are 
removed. The first condition (4) implies 2N linear equations, without Rq and T§. The second 
condition (4) implies 2N + 1 nonlinear equations, including 7}? — Rq + vo/ ft). Finally, we get 
a (4N + 1) x (4N + 1) nonlinear system with first-order or second-order polynomial entries 

F(X)=0, X = (Y,R a v T l a ,R b ,T l b ,...,R b N ,T») T , F=r fl -^ + ^. (14) 

To solve (14), three cases may be considered: 

• linear regime, for all N: (14) becomes a linear system whose solution is easy to compute 
analytically. In this limit-case, Y = and R„'° = T„' D = if n > 2; 

• nonlinear case, N = 1: (14) can be solved analytically. A detailed study shows that the 
solution may not be unique: if vo is lower than a critical value % there exists only one 
real root; if vo > % tnere are three real roots. Moreover, Y > 0; an approximation of this 
jump recovers the value Y\ deduced from the PA (11); 

• nonlinear case, N > 1: (14) is solved numerically by the Newton-Raphson method. The 
determination of F and of its jacobian J has been automatized with computer algebra 
tools: computing F and J with N = 100 roughly takes 1 second on a Pentium 3 GHz. 
The root of (14) may be not unique, like with N = 1 and vo > vo- The initialization must 
therefore be done carefully, e.g. using the exact values of the first 5 components of X at 
N = 1, and setting the other components to zero. This simple initialization works up to 
vo ~ vo- With stronger nonlinearities, this approach is coupled with a basic continuation. 

3.2 Network of cracks 

Analytical approach. A Bloch-Floquet's analysis [8] is applied to an infinite and periodic 
network in linear regime. With constant parameters, the dispersion relation is 

pea 

coskh = coskh sinkh, (15) 
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where k is the effective wavenumber and h is the spacing between cracks. If 

» = (^f, 

the waves are not attenuated. Otherwise, the waves are evanescent with decay 

3mfc = — i cosh -1 (coskh — sinttj . (17) 

Numerical approach. The analysis is much harder in nonlinear regime, since successive har- 
monics generated across cracks may not belong to the same pass-band structure. Consequently, 
nonlinear regime as well as non-periodic configurations are investigated by direct numerical 
simulations (DNS). A fourth-order ADER scheme solves (1) in the time domain. The jump 
conditions (4) are enforced numerically by an interface method [7]. With high nonlinearities, 
space-time mesh refinement is implemented around the cracks to discretize correctly the stiff 
fronts. Obviously, this approach also works with a single crack, hence it will be used to get 
reference solutions in each experiment. 



1-0 



< cos kh < 



i_-e 
T+e 



4 Numerical experiments 

A crack at a = 0.5 m in aluminium is studied: po = Pi = 2600 kg/m 3 , cq = c\ = 6400 m/s, 
K = 10 13 Pa/m, 5 = 3 10~ 7 m. Three amplitudes vo are considered: 0.01 m/s, 0.13 m/s, and 
0.6 m/s. The source is at x s = 0.25 m, and / = 100 kHz. 

First, the elastic stress transmitted across the crack is shown in figure 3. In the left col- 
umn, CT7 computed by HBM (o) is compared with the DNS (-). In the right column, the nor- 
malized harmonics 17 ((' = 1, 5) are shown. The amplitude vo = 0.01 m/s (a-b) is too small 
to mobilize the nonlinearity of the crack, and N = 2 harmonics are sufficient. We measure 
I2 = 0.017, to compare with 75 = 0.016: the approximation (12) is good. With vo = 0.13 m/s 
(c-d), N = 10 harmonics are used, and = 0.24 is measured, to compare with 72 = 0.25. 
Lastly, with vo = 0.6 m/s, N = 50 harmonics are required (e). We measure = 0.60, to com- 
pare with 72 = 1-02: the approximation (12) becomes poor. In the three cases and as deduced 
from the PA (section 3.1), the nonlinearity increases with vo- At a given vo, the I/'s decrease 
strictly with ('. 

Second, the influence of vo on the jump Y is illustrated in figure 4. The left column shows 
snapshots of u computed by DNS. The jump between mean values of u yields Y . In the right 
column, the time history of Y(t) = [u(a, t)] is computed by fourth-order Runge-Kutta in- 
tegration of (5) (-) and by PA (o). At small t, PA differs from RK 4 because it does not 
compute transients. With vo = 0.01 m/s, Y = 5.42 10 m/s (to compare with the approxima- 
tion F] = 5.026 10 -11 m/s) is not visible (a), and N = 2 in the PA (b). With v = 0.13 m/s, 
Y = 5.84 10~ 8 m/s is measured (c), to compare with Y\ = 5.60 10~ 8 m/s, and N = 6 is required 
(d). In that case, the maximum value of Y(t) is roughly 0.995: with greater values of vo, then 
max \Y\/S > 1, and the PA may not converge. It is what happens with vo = 0.6 m/s: only RK 
4 is shown (f), and Y = 9.82 10~ 7 m/s is measured, to be compared with Y\ = 1.02 10~ 7 m/s. 

Parametric studies performed by the HBM (N = 50) are proposed in figure 5. A log-log 
scale shows that Y is very close to the line with slope 2 deduced from (11), even at high vo (a). 
The amplitude of 17 increases strictly with vo, and tends towards 72 at small vo (b). 

Lastly, propagation across a periodic and finite network of cracks is simulated in figure 
6. Two spacings are considered: h = 1.72 10~ 2 m and h = 5.17 10~ 2 m. In linear regime and 




Fig. 3. Transmitted stress, with vo = O.Ol m/s (a-b), vo = 0.13 m/s (c-d), vo = 0.6 m/s (e-f). 
Left column: values of Or on a period; right column: normalized harmonics I] (i = 1, 5). 
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Fig. 4. Left column: snapshots of u around a (the dotted horizontal lines denote the mean 
value of u, on both sides of a). Right column: time history of Y = [u(cc, t)] (the horizontal 
dotted lines denote ±8). (a-b): vq = 0.01 m/s, (c-d): vq = 0.13 m/s, (e-f): vq = 0.6 m/s. 




Fig. 5. Parametric studies in terms of vq: Y (a), normalized harmonics i7 (i = 2, 5) (b). 



with an infinite network, the Bloch-Floquet's analysis predicts respectively a pass-band and a 
stop-band behavior, observed with small amplitudes (a-b). The attenuation measured in (a) is 
in good agreement with the theoretical attenuation (17). With higher amplitudes (c to f), the 
behaviors are maintained. 



5 Conclusion 

The main results of this work are as follows: 

1. with one crack, two phenomena are induced by the nonlinearity: mean dilatation of the 
crack, generation of harmonics, that both increase with A and ^ ; 

2. with many cracks, simulations show that properties of infinite linear networks are valid 
in nonlinear regime with finite networks much greater than wavelength. 

Three directions are distinguished for further investigation: 

1. investigation of shear effects, coupled or not with compressional efforts [1,9]; 

2. analysis of the periodic solution of (5), to prove Y > and dY/d > 0; 

3. effective properties of a random linear and nonlinear networks [5]. 
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